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We compute the electrical conductivity for liquid hydrogen at high pressure using quantum Monte 
Carlo. The method uses Coupled Electron-Ion Monte Carlo to generate configurations of liquid hy- 
drogen. For each configuration, correlated sampling of electrons is performed in order to calculate a 
set of lowest many-body eigenstates and current-current correlation functions of the system, which 
are summed over in the many-body Kubo formula to give AC electrical conductivity. The extrapo- 
lated DC conductivity at 3000 K for several densities shows a liquid semiconductor to liquid-metal 
transition at high pressure. Our results are in good agreement with shock-wave data. 

PACS numbers: 71.22.-fi, 72.15.Cz, 02.70.Ss 



Liquid hydrogen at high pressure has been the subject 
of intense experimental and theoretical research, because 
of its special place in the periodic table and its cosmic 
abundance. Understanding its behavior under high pres- 
sure and high temperature is important for revealing the 
properties of giant planets such as Jupiter and Saturn. 
Metallization of liquid hydrogen at high pressure is of 
particular interest. Experiments using shock waves have 
found a metallization transition'!']; at a pressure of 140 
GPa and temperature of 3000 K, liquid hydrogen has 
been reported to turn from an insulating fluid into a 
metallic one vifith a DC conductivity of about 2000 (51 
cm)~^. However, theoretically, such a metallization pro- 
cess has not been very well understood. Mean-field den- 
sity functional theory (DFT) calculations Q and Quan- 
tum Molecular Dynamics (QMD) calculations 01 have 
been used to calculate the electrical conductivity in liq- 
uid hydrogen, but these methods neglect correlations be- 
tween electrons. Such correlations are likely to be impor- 
tant for the accurate determination of transport proper- 
ties. In this letter, we propose and apply an alternative 
ab-initio method which combines the Coupled Electron- 
Ion Monte Carlo (CEIMC) method [1] with the Correla- 
tion Function Quantum Monte Carlo (CFQMC) method 

h 7]. 

Within the CEIMC approach 0, Q , the electrons and 
protons are simulated with separate but coupled Monte 
Carlo simulations, taking advantage of the separation 
of time scales in the Born-Oppenheimer approximation. 
The Born-Oppenheimer energy surface for the protons 
is determined using reptation quantum Monte Carlo 
(RQMC) |9|]. After the proton system reaches equilib- 
rium, we record uncorrelated samples of proton config- 
urations, which are subsequently used to determine the 
electrical conductivity. 

The many-body Kubo formula [13] for electrical con- 



ductivity is 
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where a = x,y, z; e (m) is electron charge (mass), V, is 
the volume, /3 is the inverse temperature, uj is the fre- 
quency, and is the a component of the momentum 
operator for the ith electron. \n) and £"„ are many-body 
eigen-states and eigen-energies of the Hamiltonian, and 
Z = Q-^En ^Yie partition function. 

The key challenge in evaluating the Kubo formula (Eq. 
[T]) is to determine the sum over all many-body eigenstates 
of the system. We compute a number of the lowest- 
energy eigenstates, and the relevant matrix elements, us- 
ing the CFQMC method Q , which we now explain. Con- 
sider the subspace spanned by a set of M many-body 
basis states fj{R), i.e., <&i(R) = J2jLidijfj(R), where 
R represents electronic coordinates. The ground basis 
state, /o(R), has Slater determinants of Kohn-Sham sin- 
gle electron orbitals times a Jastrow pair correlation and 
backfiow corrections!^ [l3]. For the excited basis states 
/j(R),j > 0, we use low-lying particle-hole excitations 
with respect to the ground state /o(R). 

Within this subspace, upper bounds to the exact eigen- 
values can be found, by minimizing the Rayleigh quotient 
with respect to c?.y , 
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where H is the electronic Hamiltonian. Furthermore, an 
improved basis set {/^(R)} can be obtained by applying 
an imaginary-time projection to the basis set: 
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FIG. 1: (Color online). Electrical conductivity (a), the cor- 
responding sum rule (b), and five lowest energy levels (c) as 
a function of imaginary-time projection t for — 1.40. Note 
the decrease of ct{lj) at low frequency is an artifact due to the 
gap in a finite system. The simulation has M = 24 excited 
states. The results are averaged over 108 twist angles (see Fig. 
131) and 10 decorrelated proton configurations drawn from the 
thermal Boltzmann distribution sampled by our CEIMC cal- 
culations. Lines are just guides to the eye. 




FIG. 2: (Color online). Electrical conductivity (a), the cor- 
responding sum rule (b), and first three lowest energy levels 
with energy difference Ei — Eo in the inset as a function of 
M, the number of trial wavefunctions in the basis set. Here 
t = 1.02 Ha"^ for = 1.40. 



Note that we have used Eq. ^ and inserted complete 
sets of basis states {|R)}- These matrix elements are 
calculated all at once with RQMC. 

Solving Eq. Q yields the many-body eigenvalues and 
eigenfunctions that are required for the Kubo formula. 
The current-current correlation functions are computed 
at half of the projection time t/2 in RQMC Q. Fi- 
nally, the S function in the Kubo formula is broadened 
by a Gaussian function during the numerical calculation, 
whose width is of the same order as the observed many- 
body energy spacing. 

We use the ground state as guiding wavefunction dur- 
ing RQMC simulations 17[ . The method is able to deter- 
mine the lowest-lying states of the system, typically fewer 
than 50 in our simulations, after which the calculation 
becomes less efficient due to the increased fluctuations in 
the matrix elements in Eqs. (5-6), a problem very much 
related to the well-known fermion sign problem of QMC. 
An analysis showsQ that the Monte Carlo (MC) effi- 
ciency must decrease with projection time as exp(— at), 
making convergence in t problematical: we call this the 
"efficiency problem" . However, it is important to note 
that as the basis is increased in size, the projected low 
energy states become more accurate Q. This suggests 
that we should include as many excitations as long as 
the efficiency is not severely reduced. 




FIG. 3: (Color online). Electrical conductivity calculated at 
the r point compared with the twist-averaged one. 



where t is the projection time. Replacing {/j(R)} with 
{/j(R)} and minimizing with respect to dij, one ob- 
tains the many-body generalized eigenvalue equation: 

M M 

^^^At)dk,{t) = Ak{t) s,,{t)dk,{t), (4) 

j=i i=i 

where Ak{t) is the fcth eigenvalue, and the Hamiltonian 
H and the overlap matrices S are defined with respect to 
the projected basis set {/^(R)} as 

U,, = J dRidR2/,(R2)(R2|ffe-*''|Ri)/.(Ri), (5) 

Sy = / dRidR2/,(R2)(R2|e-*^|Ri)/.(Ri). (6) 



All simulations reported here were at 3000 K, the tem- 
perature of the shock experiments measuring conductiv- 
ity m . We first investigate the role of RQMC imaginary- 
time projection t. Fig. [T] (a) shows the electrical con- 
ductivity of liquid hydrogen as a function of t. Defining 
S{u}) = Jq a(u) )duj , the conductivity sum rule is 
liniij^oo <S'(w) = 1 [101 • Here rie is electron density. In 
the present method, we cannot include all many-body 
states, so we do not expect the sum rule to be satisfied. 
However, it provides an indication of calculation qual- 
ity. The main measure of appropriate projection time is 
the convergence of low- frequency cr(a;). We can see from 
Fig. [T] (a) that low-frequency cr(w) curves converge be- 
tween the interval t G [0.82, 1.02]Ha~^. Here we choose 
the upper bound of t 1.02Ha^^, since we want to have 
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FIG. 4: (Color online). Electrical conductivity (scattered 
points) along with the Drude fits (solid lines) (upper panel) 
and DOS (lower panel) as a function of density for liquid hy- 
drogen. The red dot at a; = indicates the experimental DC 
conductivity value after the liquid hydrogen is metallized by 
high pressure at 3000 K [2| • The drop of electrical conductiv- 
ity to zero at zero frequency for = 1.40, 1.50 and 1.60 is 
due to finite-size effects. 



as large t as possible before the efBciency is reduced to 
an intolerable level. Note that the decrease in a{ijj) at 
w < 2 eV is an artifact due to the finite number of atoms 
in the simulation, as we discuss below. Fig. [T] (b) also 
shows that the best sum rule satisfaction is achieved at 
t = 1.02Ha-\ where S{uj) - 0.6 as — J> oo. As a fur- 
ther test of the t value, we show in Fig. [T] (c) five lowest 
energy states as a function of projection time. We see 
that indeed after t = 1.02Ha^^, the statistical variances 
increase substantially, and the estimates of the lowest 
energies are pushed down due to the noise 0] . 

Fig. [2] shows (a) the electrical conductivity, (b) the 
corresponding S{uj), and (c) the first 3 energy levels at 
fixed projection time as a function of M, the size of the 
basis set. As before, we choose M as large as possible 
while still obtaining a reasonable MC efficiency. To check 
the efficiency, we collect the same amounts of stochastic 
data for the Hamiltonian and overlap matrices in Eq. [5] 
andini solve the generalized eigenvalue equation Eq. SI 
and look at the energy gap between ground state Eq and 
the first excited state Ei . A sudden increase of the energy 
gap indicates the reduction of the efficiency, since noise 
can push down the most the lowest energy level • The 
inset of Fig. [5] (c) shows an increased slope for Ei — Eq 
curve at M — 36, indicating such a reduction of efficiency 



The corresponding a{u}) curve is also suppressed in the 
low frequency region. See Fig. [2] (a) . We find that M = 
24 is a good compromise between accuracy and efficiency. 
We also notice from Fig. [2](b) that S{lu) — 0.6 as w oo 
for M = 24; while for M = 12 and 36, the corresponding 
sum approaches 0.5. 

We note that our simulations are subject to finite-size 
errors in evaluating thermodynamic properties. Such ef- 
fects can be reduced by using twist-averaged boundary 
conditions In Fig. [3] we compare (j{uj) calculated 



at the r point of the Brillouin zone (BZ) to the one av- 
eraged over a 6 X 6 X 6 grid in the twist angle space 
(the BZ). We find that F point sampling overestimates 
the conductivity value. However, finite-size effects have 
not been entirely eliminated and are responsible for the 
observed vanishing of cr(w) at small w for all systems 
studied. This will not happen for a sufficiently large 
system in the metallic phase (e.g. at < 1.50). A 
common procedure[3, [13, 113] to estimate the DC con- 
ductivity is to discard a few data points and extrapolate 
the higher frequency data to lu — . For systems in- 
side the metallic phase one can use the Drude formula 
a{uj) = (Tdc/(1 + w^r^) to estimate the DC conductivity 
ctdc Sind the electronic relaxation time r. 




FIG. 5: (Color online). Electrical conductivity (scattered 
points) and Drude fits (solid lines) as a function of density 
for liquid hydrogen. The red dot at a; = denotes the ex- 
perimental DC conductivity value after the liquid hydrogen 
is metallized by high pressure [ij]. 



We now proceed to study the electrical conductivity 
at various densities at T= 3000 K, where the metalliza- 
tion of liquid hydrogen has been observed [ij. We run 
CEIMC simulations at r, = 1.40, 1.50, 1.60, and 1.65 
for 32 protons and electrons. Results (scattered points) 
along with the corresponding Drude fits (solid lines) are 
shown on the upper panel of Fig. |3]and many-body den- 
sity of states (DOS) on the lower panel, where all the 
energy values are calculated with respect to the ground- 
state energy of the configuration. We can see that DC 
conductivity decreases from about 8600 {il cm)~^ (see 
below) at = 1.40 to about 3500 {fl cm)^^ at = 1.60 
and to nearly zero at = 1.65. The dot at uj — 
in Fig. 0] shows the experimental DC conductivity after 
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metallization The conductivity points for Ts = 1.65 
near zero frequency has typical features for liquid semi- 
conductors, while a 2000 (fi cm)^^ is a typical liquid 
metal value [l5| . Thus, based on 32-atom simulations we 
see a liquid metal to a liquid semiconductor transition 
in a density region of 1.60 < < 1.65 at T=3000 K 
for liquid hydrogen under high pressure, which is close to 
the experimental transition density of = 1.62. Such a 
metal-semiconductor transition is also hinted in the DOS 
figure. Going from — 1.40 to = 1.65, one can see 
the gradual opening of an energy gap. 

We further look at the metallic behavior of the liquid 
hydrogen at = 1.40. The Drude fit gives ctdc = 8620± 
1000 {VL cm)-i andr = 2.0±0.2a.u. or 3.1±0.3x 10"^^. 
A rough estimate of electron mean free path {I) then gives 
I ~ 1.2 A, which is about 1.5 times the average proton- 
proton distance (1.4 Bohr) necessary for the system 
being metallic. For — 1.50 and 1.60 we also fit o'(w) 
points to the Drude formula. The fitted curves become 
flatter as the density decreases, which signals that at Ts — 
1.60, liquid hydrogen has a small electronic relaxation 
time, making it a bad metal, similar to liquid carbon 

Finally, we explicitly test the importance of finite-size 
effects by performing simulations for 54-atom cells as 
shown in Fig. [5l The energy gap at = 1.65 that is 
present in the 32 atom simulation now disappears. The 
liquid semiconductor [at = 1.75 the extrapolated DC 
conductivity is around 1300 {Q, cm)~^ less than a typi- 
cal liquid metal value 3, hence the name] to liquid 
metal transition density is estimated to be in the range 
1.65 < Ts < 1.75, which is again close to the experimen- 
tal density of = 1.62 IJ. Further increase of system 
size is not possible at present, but previous CEIMC cal- 
culations have found that the ground-state energy with 
54 and 108 atoms are very close if twist-averaged bound- 
ary conditions are used. Therefore, we expect that the 
liquid metal transition density determined above is close 
to the thermodynamic limit. A more definitive answer 
will require significant additional calculations. 

We have proposed a method to calculate the frequency 
dependent electrical conductivity using CEIMC with cor- 
related Monte Carlo sampling and the many-body Kubo 
formula. The method is able to estimate some of the 
lowest-lying energies and their corresponding overlap ma- 
trices, and hence it is suitable for the extrapolation of 
DC conductivity from the AC conductivity. With this 
method we study the DC conductivity of liquid hydrogen 
as a function of density at 3000 K and at high pressure, 
and show the metallization. Our DC conductivity values 
at the metallization point is in good agreement with the 
shock- wave experiments [ll] . Furthermore, our metalliza- 
tion density (r^ 1.65) is close to the experimental value 
of Ts — 1.62. In the future, it will be interesting to apply 
the method to more densities and temperatures to have 



a complete understanding of liquid hydrogen electronic 
transport properties. It is also very promising to extend 
our method to other systems, such as lattice models. 
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